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Abstract 

Nonlinear transport through a quantum dot is studied in the limit of weak and strong intra-dot 
Coulomb interaction. For the latter regime the nonequilibrium self-consistent mean field equations 
for energies and spectral weights of one-electron transitions are formulated. The increase of the 
bias-voltage window leads to a strong deviation from Gibbs statistics: the populations of states 
involved into a tunnelling are equalizing in this limit even at low temperature. For a symmetric 
coupling of a quantum dot to two leads we provide simple analytical relations between heights of the 
current steps and degeneracy of a spectrum in a two-dimensional parabolic dot in a perpendicular 
magnetic field in the both regimes. 

PACS numbers: 73.23.-b,73.21.La,75.75.+a 
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I. INTRODUCTION 



Nowadays, quantum dots (QDs) are among common experimental tools to study proper- 
ties of correlated electron systems at atomic scale, in a controllable manner. In particular, 
transport measurements through QDs provide a rich information on internal dynamics gov- 
erned by external gate and bias voltages, and by the degree of coupling of the dot to the 
source and drain electrodes. Depending on the above conditions, a spectrum of a quantum 
dot, extracted from the single-electron spectroscopy, manifests a variety of features from 
shell structure to chaotic dynamics (see for review RefsQ, I2). 

One of the simplest quantum devices to study transport consists of two metallic contacts, 
attached via insulating/ vacuum barriers to a QD. Combination of pieces with very different 
physical properties suggests that theoretical descriptions may require different approxima- 
tions in accounting for the interactions in the constituents. In addition, the openings affects 
an intrinsic structure of initially closed quantum dot. Evidently, the energies of many- 
electron transitions (levels) acquire a width and both single- and many-particle processes 
contribute to the life-time of the levels. Moreover, the levels can be shifted due to the 
coupling of the QD to the leads. The shift can also be spin-sensitive due to many-electron 
nature of the electron transitions involved jjj. 

For small QDs, when the size quantization of the carrier motion becomes important, 
and at low enough carrier density the strong electron correlations (SEC) affect the electron 
transport through the dot. At the weak dot-lead coupling and relatively small temperatures 
(< 200 niK), the intra-dot SEC are dominant in the tunnelling. At this (Coulomb blockade) 
regime, the charging energy E@, that needed to add an electron to the dot with N electrons, 
considerably exceeds the level spacing and the thermal energy ksT. According to the clas- 
sical Coulomb blockade picture, the transport is only allowed when the N and (N+l) states 
are both energetically accessible (cf Ref|J). The conductance usually shows a Breit-Wigner 
resonance as a function of gate voltage indicating that electrons pass by a discrete level in 
the quantum dot by a resonant tunnelling. At lower temperature the co-tunnelling events 
predicted in Ref 5] become visible; they involve the simultaneous tunnelling of two or more 
electrons over virtual states, giving rise to the current inside the Coulomb diamond 
With a further decrease of the temperature, Kondo-type effects can develop if the system 
has quasi-degenerate localized states. For example, a weak current through QD, observed 
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a, 



within the Coulomb blockade regime |7j], is interpreted as a Kondo effect (see for a review 
RefQ). 

Recently, in experiments with a small quantum dot in the Coulomb-blockade regime 
a fine structure was observed in the conductance as a function of gate voltage versus source- 
drain voltage. It was suggested that this phenomenon is mainly due to a co-tunnelling jsj]. 
The theory of co-tunnelling p neglects, however, specific quantum effects caused by a strong 
Coulomb interaction (SCI), which affect the transition energies and the tunnelling rates. On 
the other hand, numerous papers devoted to quantum effects in transport through QDs are 
focused on the analysis of Kondo phenomenon in single-level QD models (cf Refsllfll) in a 
linear-response regime, which is well understood nowadays. However, the single-electron 
spectroscopy clearly indicates that a shell structure of small dots plays an essential role in 
the transport and weak dot-lead coupling 4] when the temperature is still low, but above 
the Kondo temperature. This becomes especially evident, when the confining energy exceeds 
the charging energy (hereafter, this regime is called a weak Coulomb interaction (WCI) 
regime). One of our goals is to present a self-consistent approach for a nonlinear transport 
through a multilevel quantum dot in the SCI regime. We will demonstrate that even in 
this regime a shell structure of a small dot can be extracted from the analysis of nonlinear 
current as a function of a finite bias voltage. We will show also that the conductance inside 
of the " Coulomb diamond" is governed by spectral weights of one-electron transitions, that 
complement the co-tunnelling picture. 

In this paper we consider the transport through a small QD with a parabolic effective 
confining potential. As an input, we use the intrinsic structure of a closed dot. The use of 
the parabolic confinement is well justified by experimental observations, for example, magic 
numbers of a two-dimensional harmonic oscillator |l2| and a fulfilment of the Kohn theorem 
[13I ] in far- infrared spectroscopy experiments 3]. The eigen problem of interacting electrons 
in a closed parabolic QD can be solved approximately within an unrestricted Hartree-Fock 
approximation (cf RefllfJ) or local density approximation to a density functional theory 
[2]. When the confinement dominates the Coulomb interaction, the eigen spectrum is well 



approximated by the Fock-Darwin levels 



161 ] . Recent experiments nicely confirm this fact, 



indeed (see, for example, experimental results in Refllll). This allows to describe a QD 
within a single-particle picture of electrons in an effective parabolic potential (hereafter we 
denote such states as I7) ). We consider the SCI regime in the limit of very strong intra-dot 
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interaction ("Hubbard" U — > oo). In this case the Hamiltonian of QD at the SCI regime is 
diagonalized by the states |0) and (7), and only the transitions between empty and single- 
electron occupied QD states contribute into transport. The assumption that the states |0) 
and I7) are well separated by a large energy gap from two- , three-, etc correlated states 
can be easily verified, for example, within an analytical model for two-electron parabolic 
quantum dot At zero coupling to leads one-electron energies coincide in both the WCI 
and SCI regimes. It presents a nice opportunity to study in a transparent way the role of 
the Coulomb interaction within this particular charge sector. Although the bare energies 
of QDs coinside n these two regimes, in the SCI limit a hopping of an electron into a QD 
is determined by the occupation numbers of the QD states. We employ here the Hubbard 
operator technique [18] which is specially designed for such cases. The mean field theory 
in terms of the Hubbard operators is quite different from the one constructed for the WCI 
regime. For example, the consideration of a single spin flip in the infinite U Hubbard model 
shows that the mean field theory of this type gives the same resu ts with both Gutzwiller 
wave function and the three-body Faddeev equations in this limit 

As pointed out above, we will analyze the transport through QD in non-linear with respect 
to the bias voltage regime at the weak dot-lead coupling. In other words, the transparency 
of junctions is assumed to be small. The latter requirement has a pragmatic basis as well: 
enormous efforts for a fabrication of QDs are motivated by the desire to understand the 
inner structure and to exploit it in various applications. 

Thus, the coupling to the leads will be treated perturbatively. Below we focus on the 
transport in the resonant-tunnelling regime, where the Kondo physics is not involved. In the 
language of the Anderson impurity theory, we will be interested in the intermediate- valence 
regime, when QD levels are located within the range of the resonant tunnelling. We recall 
that in this region such a widely used method as a mean-field approximation in slave-boson 
theory is not applicable, in contrast to our approach. 

The paper is organized as follows: in Section II we will discuss our model Hamiltonian, 
the nonlinear current in terms of the Hubbard operators and formulate the mean field 
approximation. Section III is devoted to the derivation of equations for population numbers 
and energy shifts in the diagonal approximation. The results for the transport through a 
parabolic quantum dot in a perpendicular magnetic field at weak (WCI) and strong (SCI) 
intra-dot Coulomb interaction are given in Section IV. The conclusions are finally drawn in 
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Section V. Two Appendices contain some technical details. Some preliminary results of our 
analysis are presented in Refl^ll 



II. FORMALISM 

A. Model Hamiltonian and current 

As discussed in the Introduction, we consider the system which is decribed by the follow- 
ing Hamiltonian 

H = H l + H r + Hqd + H t , (1) 

where 

H X = £\kac\ ka C Xka 

Hqd = ^^e 7 n 7 + £7" 77 /n 7 n 7 /, n 7 = d^d^ 

7 77^7' 
ka£X,X,f 

The term H\ describes noninter acting electrons with the energy e\ka, wave number k and 
spin cr in the lead X = l,r. The closed dot is modelled by Hq D . Here e 7 is the energy of 
electron with a set of single-electron quantum numbers 7 in the confining potential of the 
closed QD at U — 0. Tunnelling between the dot and leads is described by the term H t ; the 
matrix elements v kaQry here couple the left and the right contacts with a QD. The ratio of 
level spacing of noninteracting electrons in a QD to "Hubbard U" , \e 1 — £y| fU = Seyf/U 
determines whether the intra-dot electrons are in the regime of weak (Se^y/U > 1) or strong 
(dSyf/U < 1) Coulomb interaction. 

The Hamiltonian of QD is transformed into a diagonal form by introducing the Hubbard 
operators j3] (see also Appendix EJ. As a result, the Hamiltonian ([TJ acquires the form: 



H - J2kaX £ \kac[ ka C X ka + £oZ°° + ^ E^Z 

+ Ek„xrf (<Ak* X01 + h - c ) 



(2) 
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Following Refl23l. we obtain (see Appendix EI) for the "left" steady current 

J > = % E/^{ K, 7 M-r;, 7 M]G< 7 ,H 
77' ^ 

x [G^,(u) - G^,(u)]\ (3) 

The current (jHJ) is expressed in terms of Green functions (GFs) defined for the Hubbard 
operators : 



Gl Y (co) = J dit-t'y^Gl^t'), V =<,>,R,A 
G< y (t,t')=t(x^(t')X^(t) 



G> 7 ,(M0 = -*(^ O7 (^V) 

G^(t,t>) = ^ie(±t T t') ({x^(t),x-y'°(t')}) (4) 

The GFs depend on difference of times (or one frequency u) because we consider the steady 
states only. Note that for a wide conduction band, i.e., when the bandwidth is much larger 
than any other parameter in the problem, the width function 

r 7Y H = 7T Y1 v YmK" - txko)vl n (5) 

ka 

weakly depends on uj and will be replaced by a constant when we will consider a particular 
model (see also Appendices IAIB|) 

Thus, we have to find the GFs G 7 ,. As well-known, the Dyson equation does not exist for 
any many-electron GF (particularly, for the GF defined on the Hubbard operators), since the 
(anti) commutation of many-electron operators (see Eq. (jA.12jO generates not a "c" -number 
as in case of fermions or bosons, but an operator again. The perturbation theory for the 
many-electron GFs reflects this fact by generating the additional graphs compared to the 
conventional theory for fermions or bosons (examples are given in Refill also see a Sec. lib). 



Therefore, those of results of Refs 



23 



2J, which are derived on the basis of a Dyson equation 



for the retarded and advanced GFs, should be re-inspected. As will be seen below, the 



formulas analogous to the ones derived in Refs 



23 



2J are valid only within the well-known 



Hubbard-I approximation and mean-field approximations (MFA) [18| (Refll8lwill be referred 
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below as I). The same approximation is sufficient for calculation of the transport properties 
at small transparencies of the junctions too. 

Let us introduce the full width T and dimensionless partial widths a: 

l-\'7 = + Fy , ^7'7 = ^-y^/^Vt' (^) 

When the coupling to leads is so small that the widths T of the transitions are much smaller 
than the temperature T, (3Y <C 1 (where (3 = 1/T) and the level spacing, r/5e 7 y <C 1, 
the mixing of levels due to non-diagonal terms does not influence the physics of transport. 
Namely this range of parameters is exploited usually in experiments and will be in the 
scope of our interests (starting from Sec. III). Neglecting the non-diagonal terms leads to the 

7 ; 7 ^7' 7-^75 ^77' 5y><yOlj. 

becomes: 



diagonal approximation, where Ty 7 = <5y 7 r 7 , a* , = 5 7 ' 7 a^. The expression for the current 



' Y 7» 



+ [a\m - aj;/ P H] [G*„(u) - G*»] \ (7) 



The advantage of the Wingreen and Meir formulation [2^J lies in the fact that the current 
(in our case Eq.fjHJ)) is expressed in terms of the GFs of a QD only. We recall that this result 
can be obtained only for non-interacting conduction electrons. Thus, our next-step goal is 
to derive equations for GFs G^^ A in the MFA. 

B. Exact equations on imaginary-time axis 

Our derivation is based on the approach 0, [l^ that emp 
Schwinger and developed further by Kadanoff and Baym 
the intra-dot Coulomb interaction, the conduction-electron subsystems are linear and can be 
integrated out. Instead an additional effective interactions between electrons in a QD arises. 
Following RefQ, equations for the QD subsystem will be written in terms of functional 
derivatives with respect to external auxiliary sources. The derivation of the equations will 
be performed in three steps. Firstly, we derive exact equations for the QD GFs on the 
imaginary-time axis. Secondly, the MFA will be formulated. Thirdly, we will perform a 
continuation of these equations onto the real-time axis. 



oy s the ideas suggested by 
26|. Since we consider only 



The time evolution of the Hubbard operator X° 7 is defined by the commutator 

[x°\h] = A^x°-y+ (8) 

+ E <M 0W°° + Z^) c Xk „, 

which contains non-linear terms. These terms produce GFs of the type 
-i(T (5 yi Z 00 {t) + Z^{t)) 

c \k"a(t) c xk'a(t'))- Below we omit the summation signs: 
summation over repeating indices is implied unless otherwise stated. 

To proceed further, we introduce the auxiliary sources U 00 (t), IT™ (t) into a definition of 
the GFs: 

i u i {TS{t ,p)) 

'u 00 (t 1 )Z 00 (t 1 ) + 

(9) 



(rto-ip 
— i I dti 

+J2u 77 '{t 1 )z' r ' r '(t 1 



77' 



Here A,B = c\ka, X° 7 ; (A = I,) for for the left and (A = r) for the right contacts. We will 
use the following notations: 



<Wf,tf) G^,(t,t>) ) * \ (TX^(t)c{, k>tT ,(t')) u (TX^(t)X^(t>)) 
The equation of motion for the GFs Ck a ,k'a'(t,t') and Gfe CT)7 '(t, if) are simple: 



(10) 



+ ^< y ,GV, fc v(M') (11) 



(id t — S\ ka ) G ka) y(t,t') — V k X arfl >G-y",'y'{t, t') 



(12) 



As seen from Eq. fjTTj) . a zero GF C^ fc/(7 , satisfies the equation 

/ dhY.m-exkMt-hnc^A^t') = (is) 



$kk'8<r<r'8(t - t') 



Hereafter the integrals are taken between to and to — iP, if the other limits are not specified 
explicitly; here to is an arbitrary moment. With the aid of this equation, Eqs.([llj ) .()12 p can 
be cast in the form of integral equations: 

Ckcr,k'cr'(t,t') = C^J(t, t')8 kk l5 a(T ' + (14) 

72 

(15) 

Analogously, using the equation for the GF Gy^v^ijO with respect to the right time t', 
we obtain 

GV,*v(M') = E / ^iG 7 » )7l (Mi)<iv^(ii,0 (16) 

71 J 

The equation for G^(t,t') = -i {TX^{t)X^' Q {t')) u has, however, more complex form: 
[S^idt - Aj(t)] G lia ,{t,t') = 8{t - t')pn\t) + 



fc(rGA,A,7i 



(17) 



where 



P^\t) = ^T|x° 7 (t),X yo (t)}^ (18) 
AW ( t ) = Af5 717 + [Urn® _ 5lll U°°(t)] 

Here A? = Ey — So is the transition energy between the single-electron state and the one 
without electrons in a QD. We choose e = 0. The term (U^' (t) - <5y 7 £/ 00 (t)) is due to the 
differentiation of S(t , P) (cf Ref 26). 

Below, we will proceed under the standard assumption that the interaction is absent at 
infinitely remote time to — > — oo and is switched on adiabatically. This assumption does not 
affect our approach, since we study a stationary regime. In this case we use the relation: 

(TZt(t)X^(t)B(t'))u = 
((TZ*(t)) + ^smm) (TX*y(t)B(t>))v (19) 
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that follows from a functional derivative 5/5U^(t) of Eq.® (here £ = [00], [77']). 

With the aid of Eq. (jl9|) . using Eq.([15j) and the definition for the effective interaction via 
conduction electrons 



Here, we used also that SCj^J/SU = 0. 

Eq. (j21J) is the exact equation for the GF G ia >(t,t') on the imaginary time axis. Since all 
other GFs are expressed in terms of the G 7 », 7 , the iteration of Eq. (j21)) with respect to the 
effective interaction generates a full perturbation theory. However, within this formulation 
of the theory, the continuation from the imaginary time axis to the real one should be 
performed in each term of the expansion. 

The well-known "Hubbard-I" (HI) approximation can be obtained from Eq. (}2*Tj) by 
putting therein SG^i-JdU = 0. Let us define a zero locator D^,(t,t') 




(20) 



(see also Eq. (jB.7jO . we transform Eq. (fT7j) to the form 



[5, yi id t - A% (t)] G 71 ,y(M') = 6(t - t')P^\t) 




(21) 




5(t- t')5, 



(22) 



In the HI approximation, the locator is defined by a solution to the equation: 




(23) 



where 




(24) 



7l> 
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Using the above equations, we rewrite Eq.(|21|) in the form 



-£ 7)71 (Mi) Ig 71iY (^,0 = 6{t - t')P^\t') 



with 



(25) 



/ dtiE 71 S 7 ™(t,*i)G w (ti,t') = (26) 



One may notice that Eq. (|2fij) has a simple form in symbolic notations 

D H )G + V^ = P (27) 
This equation suggests that the full locator D can be defined by the equation 

£ / dhlD^t^G^yi^t') = 5(t - t')P^'(t'), (28) 

71 J 

where 

and the full self-operator is defined as 

5* 7i7l = S^ ni + S^ r/1 (29) 

The magnitude 5 7)7l has been named in I as "a self-operator" in order to distinguish it 
from the standard self-energy operator. In fact, the magnitude S^ trL (t,tf) renormalizes the 
transition energies in a QD as well as the end-factors P 77 '(t). Indeed, one can see from 
Eq. (j2B|) that 

G 71>Y (t,t') = D^ 72 (t,t')P^(t') (30) 

and, therefore, P 727 '(£') enters into the equations for the GFs under the sign of the functional 
derivative. Being an expectation value from the anticommutator of many- electron operators, 
it is not a constant (like in the cases of Bose- or Fermi-operators). Therefore, its functional 
derivative is nonzero and generates a sub-set of graphs that do not appear in standard 
techniques for fermions and bosons. These graphs describe kinematic interactions. 
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C. The energy shifts (imaginary time). 



An effective interaction of the intra-dot states via conduction electrons arises in both the 
WCI and SCI regimes. In both cases the lowest orders of perturbation theory contains the 
terms with integration over wide energy region of the continuous spectrum and, therefore, 
it results in finite widths of the intra-dot states. Particularly, in the SCI regime the width is 
provided already by the Hubbard-I approximation. Below we will show that in contrast to the 
WCI regime, in the SCI case the width depends on the population of the many-electron states 
which are involved into single-electron transition in question. This fact is more or less trivial 
consequence of non-Fermi/Bose commutation relations between Hubbard operators. A non- 
trivial consequence of it is kinematic interaction that shift transitions energies in a QD. 
The correlation-caused shift arises in a natural way within the mean-field approximation in 



slave-boson theory (see, e.g., Ref. 27]). This approach misses, however, the spin and orbital 
sensitivity of the shifts for different levels. This obstacle can be resolved by means of the 



Hubbard operators formalism (cf I and Refyl). In the present paper we use also the mean 
field approximation, which is similar in spirit to the one of Refl27l The details and physical 
meaning of our MFA, however, are different. 

The MFA for the Hubbard operator GFs have been defined in I on the imaginary-time 
axis by two conditions: first, the dynamical scattering on Bose-like excitations is neglected, 
i.e., SP/SU = and, second, full vertex T = SD^/SU is replaced by zero one, r = 
SDq 1 /5U . In other words, the MFA omits the fluctuations around the stationary solutions. 
The advantage of formulation of the theory in terms of Hubbard operators compared the 
slave-boson technique consists in the following. The MFA can be derived for a general case 
of an arbitrary Hubbard-operator algebra, or, in physical terms, for the case of arbitrary 
number and a type of (Bose- or Fermi-like) transitions (see I). Besides, one can go beyond 
the MFA in a systematic way. For example, in the presence of the long-range part of the 
Coulomb interaction one can reformulate theory in terms of screened interaction (see Ref. [3] 

)• 

In order to develop the MFA we have to extract from the self-operator S the local- 
in-time contribution in the lowest order with respect to the effective interaction V. The 
self-operator 5* (see Eq. fl26J) ) is already proportional to the interaction V 7J i. Therefore, one 
should calculate only the derivative 6G/ 5U in the lowest order with respect to the interaction. 
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We recall that the interaction Vy 7 (t, t\) does not depend on U^. Thus, the derivative is 

SG„ y(ti,t 2 ) <5[0 72 , T5 (tl,t2)PT57' '( t2 )| 



sw(t) sm(t) 
The second term in Eq. (J31j) expresses the dynamical re-scattering process 



(31) 



Spntf/su* oc (TZ 00 ZZ) - 
- (TZ 00 ) (TZ*) + <TZ 7475 Z«) - (TZ 7475 ) (TZt) , 

where the Fermi-like excitation, described by the locator D, scatters on the Bose-like 
excitation, described by the GFs 5P/5U. The diagonal correlators (TZ 77 (t)Z 77 (t')) — 
(TZ 77 (t)) (TZ 77 (f )) and (T Z 00 (t)Z 00 (t')) - (TZ°°(t)) (Z 00 (t')) describe the fluctuations of 
the population numbers and (TZ 77 ) = iV 7 and (TZ 00 ) = N , respectively. The non-diagonal 
GFs (TZ 11 ' (i)Z 1 ' 1 (t')') describe the transitions between one-electron states 7 and 7'. We 
consider the resonant tunnelling (non-Kondo regime !) at low temperature and at small 
transparencies of junctions. There the fluctuations of the population numbers produced by 
the dynamics can be ignored in this region of parameters. This makes the MFA applicable. 
Taking into account only the first term in Eq. (j31jl . we obtain 



S J}73 (t,t 3 ) =i^2 / dt i I dt 2 V J1J2 (t,t 1 ) x 

7172 



x I &rn su oo {t) + SU71j{t) | G W***0 



7172 " 



^- ^72,75(^1^ ^ 2) 
5U 00 (t) 



+ ^*^ (32) 
The derivative of the locator SD/SU can be calculated using the trick 

5D/5U = -D (5D- X /5U) D, (33) 

that follows from the equation 5 (D~ X D) /SU = 0. The validity of the trick is well justified, 
since we consider only stationary processes and within the assumption that the interaction, 
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which is absent at remote time to — > — oo, is switched on adiabatically. By means of 
Eqs.dnn},®, we have 



S ir/ '(t,t f ) ~ —i J dt\ J dt 3 V yil2 (t,ti)D^ 2 ^ 3 (ti,t 3 ) x 
x (S^SD-^t^/SU^t) +5D-] y (t 3 ,t')/5U^(t)) 



Replacing the full vertexes 

<5C/«(t) 

by zero ones and using Eqs. (fTHJ) , (J221) , we find: 



Tf' y (t 3 ,t';t) = ! (34) 



rir (ts,h;t) ~ [rW(t 3 ,ti;t)]^ 7 = 7o 7 o' 7 (*3,«i;*) 



-£(t 3 - 



<«7 00 (i) v ° ' 5U 00 (t) 

~ m ~ tl) WW) 

5 73 y5(t 3 -ti)5(t-ti); (35) 



l^'(t 3 ,t l] t) = -6(t 3 -t 1 )x 
5[A( 3 ) 5 y7 3+[f/ 737 '(ti)-5y 73 f/ 00 (ti)]] 

X 



5U^{t) 

= — 5 737l 5 77 '5(t3 — ti)S(t — ti). (36) 

As a result, we obtain 

Sytf(t,t?) = -i J c^^ 1 y 772 (^ ) ^l) J D 72I y(^ 1) ^ + )5(^-^ , ) 

-H<5 771 I dt 1 V lir2 {t 1 ti)D 12)11 {ti ) t + )5{t — t r ). (37) 



As seen, this contribution is indeed local in time and, therefore, represents an effective field, 
which shifts the transition energies in the quantum dot. 



III. THE DIAGONAL APPROXIMATION 

The analysis of magneto-transport through a multilevel quantum dot, even within the 
MFA, is still a formidable numerical task. Fortunately, as discussed above (see the arguments 
between Eqs.© and (J7J)) we can use a diagonal approximation. This makes the problem 
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much easier numerically and, what is more important, it allows to get advanced quite far 
in analytical treatment. The non-diagonal contributions certainly are important for the 
formation of Kondo state, however, their role is less important in the resonance-tunnelling 
regime. 

In the diagonal approximation the self-operator takes the form 



x(l-5 77l ) (38) 

As seen, the self-interaction is cancelled. This nice feature of the theory is provided by 
the commutation relations. The mechanism is as follows. We recall that in the equation of 
motion for the Hubbard operator, in the first order with respect to the tunnelling matrix 
elements v, the anticommutator {X° 7 , v^X j0 } c p generates the operator P° 7 = Z 00 + Z" n . In 
order to obtain the first order in the perturbation theory with respect to effective interaction 
V = v'C(t,t')v we have to calculate the second order with respect to v. This gives rise to 



the commutator 



^07 pO-y 



[X° 7 , Z 00 ] + [X° 7 , Z 77 ] = -X° 7 + X° 7 = . Thus, indeed, 
the cancellation of the self-interaction is an inner property of the theory. 

In the diagonal approximation the Hubbard-I term of the self-operator, Eq. (j24j) . is 



S%(t,t 1 ) = I™(t)V n (t,t 1 ). (39) 

Since there is only one type of coupling between transitions, [0, 7] •<=>- [7, 0], we can use the 
simplified notations: 

CHI — cHI 077 — p T r — T/ 

A 7 ° 7 } = A 7 0) , ^ Eff\ = D„ ... (40) 
Equation for the locator, Eq. (J23*j) . transforms in these notations to the form 

J dt l {d- 1 {tM)~ S^ I (t,t 1 )}D J (t 1 ,t') = 5(t-t'), (41) 

where 

d-\tM) = [idt - Af(t) - Ef f*(t)] 5(t - h). (42) 
Equation for the GF G 7 (= G 7i7 ) can be written as follows: 
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/ dh E 71 1 [tit - Af(t) - sf /*(*)] <y(t - t x ) 

-5f J (t, ti) Jo^t') =6(t-f)P^f). (43) 

Now, when the equations for the GFs and locators in the MFA are derived, we can put 
auxiliary fields U*(t) = 0. After that the Hamiltonian does not depend on time any more 
and, therefore, the locators and GFs depend only on time difference. P<y(t') does not depend 
on time as well: P 1 (t') — > P 1 = N + iV 7 . The latter simplifies the problem of derivation 
of the expressions for S^ 1 ^ — ti) for real times (R stands for Retarded and A for Advanced 
GFs): 

lS? z (t - tOf A>,< = P 7 [V^t - tl )] R > A >>>< (44) 

For the transition energy shift (that also does not depend on time at U^(t) = 0) 

the "lesser" value must be used. As a result, Eq.([41j). written in the form of the integral 
equation on the imaginary time-axis 

£) 7 (t _ t') =d 1 (t- t') + (45) 
J dU J dt 2 d 1 {t - t 2 )P 1 V y (t 2 - h) D 7 (t a - t'), 

can be immediately continued into the real-time axis. 
Correspondingly, the continued equations have the form 

D* /A {t - f ) = d* /A {t - t') + (46) 
d? A (t - t 2 )P 7 < M (t 2 - h) P? M (ii - f ) 



D<(t-t') =d<(t-t') + 

/oo poo 
dh / dt 2 [d J (t-t 2 )P~ f V J (t 2 -t 1 )} < D A (t! -t') 
oo J —oo 

/oo roc 
dh / dt 2 d*(t - t 2 )P 7 Kf (t 2 - h) D<(h - 1'), 
-oo J — oo 
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or, opening the square brackets, 

D<(t-t') =d<(t-t')+ (47) 

/OO POO 
dt x / dt 2 d<(t - t 2 )P y V A (t 2 - h) D A (tx - t') 
-co J —CO 

/CO POO 
dt x / dt 2 d R (t - t 2 )P^V<{t 2 - *0 D A {tx - t') 
-co J —oo 

/CO POO 
dt x / dt 2 d R (t - h)P 7 V R (t 2 - tO D<(tx - if). 
-oo J — oo 

Using the matrix form, we can simplify the latter equation 

[1 - d R PV R ] D< = d i? PV < B A + d<PV A D A . (48) 

As follows from Eq. (J46|) . the square bracket in the left-hand side of Eq . (I48|) is d R [D fl ] . 
Multiplying Eq . (J48j) from the left by [d^] 1 and taking into account that [d^] 1 d^ = 1, 
[d^] 1 d< = 0, we have 

[D fl ] _1 D< = PV<D A (49) 
Thus, we find the same relation which is known in the theories for fermions (see Ref 26): 

D<(t - t>) = dtx JZo dt 2 D R (t - t 2 )P 7 x 

xVffo-^Dfa-lf) (50) 

The advantage of this expression consists in the fact that it does not contain bare, non- 
renormalized magnitudes. On the other hand, its validity is restricted to: i) MFA and ii) 
stationary states (the step from the differential equation Eq.()49p to the integral one, Eq.(|5U|). 
uses the boundary condition; the same statement is valid for Eq. (14611) as well. The expression 



for D> obtained in a similar fashion, leads to Eq. (}50j) . where VJ* should be replaced by V?. 
Making Fourier transformation, we find: 

1 



D R ' A { 



d R ^ A (uj) 



-i 



cu- A^ - P 1 V R/A (u)±i5 
) 7 K ( W )P/<->(o;)^ 

Multiplying Eqs.JSU,® b Y P i we obtain GFs (see Ea. lEflty ): 



(51) 



D<>>H , D K ( W )P/<.>( W )D^( W ). (52) 



G R/A (u) = L (53) 

G<>>{«,) =G^MF<.>( W )GfH. (54) 
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The effective interaction V>' K expressed in terms of the width functions YHuj){\ = l,r; 
see also Appendix |Bj) is : 

Vy(u) = -i^^H(l-/ A (o;)) 

A 

W = iJ2 T »f^) ( 55 ) 

A 



where 



V 



r 7 (w) = r» + r», 

where P denotes the principal part of the integral. The real part A 7 (u;) in the re- 
tarded and advanced interaction Vj A (u) can be neglected in a wide-band case (A 7 (a>) 
[(|^7,fcAo-| 2 ) /W~\ In |(W + ui) I (W — uj)\, where W is of order of a half of a bandwidth and 
we are interested in uo <C W ). Then (u) ~ r RT 7 and, therefore, 

Using this equation, we obtain a few useful relations that will be used below. In particular, 
the difference of GFs is 

G*{w) - G A (cu) = -2tuP 7 L 7 H, (57) 

where 

L » - 1 *f! n (P r ,2 (58) 

(to - Ay) + (P 7 r 7 ) 

is the Lorentz distribution. Note that, contrary to the WCI case, in the SCI case the width 
of the distribution depends on the non-equilibrium population numbers via P 7 = Nq + N T 
With the aid of Eqs.fJSS},©, we obtain for the GF G<: 

G<(u) = G*(u)V<(u) G A (u) = 2™P 7 L 7 H/», (59) 

where a 7 r = r 7 r /T 7 and f(ui) is the weighted Fermi-function: 

/»= [a 7 /,(u,) + a 7 / r (u;)]. (60) 
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A. Equations for population numbers. 



All the locators and GFs depend on the single-time correlators P 1 = 
(X°yXT° + A 7 °A° 7 ) = N + JV 7 . The expectation values (X 7 °X° 7 ) and (A° 7 A 7 °) 
can be found by means of GFs G^ ,:> {t : t). We recall that our aproach is developed for the 
stationary case only Thus, we have 

<X 7 °X° 7 ) = -?(i(X 7 °(t)X° 7 (t))) = 

= -iG3t,t + ) = -i f duG<(uj) (61) 



Since the GFs depend on P 7 themselves, Eqs. (|61|) for the population numbers are, in fact, 
the equations of self-consistency: 

iV 7 =-iJ dujG<{uj) 

N =ijduG>{uj) (62) 

The normalization condition, which follows from the averaging with respect to the 
state/ensemble of interest of Eq. (jA.5|) . provides the additional equation 

l = N + J2 N i ( 63 ) 

7 

Substituting Eq.([59|) into Eq.(|52j) we obtain 

A 7 = P 7 s 7 , s 7 = J duLj(uj)f 7 (uj). (64) 

From Eqs. (jo^j) .(fo^ j) . (|6^|) one may conclude that the magnitude s 7 is the integrated "lesser" 
part of the locator: 

D<(to) = 2niL^uj)f(uj). (65) 

For Tf3 1, the Lorentzian can be approximately treated as a delta-function. Therefore, the 
integral in Eq.(jH2I) can be approximated by values of the weighted Fermi function, Eq. (|6T)|) . 
at the transition energy. It results in a simple expression: 

A 7 = P 7 S 7 , S 7 ~ [MA^ + / r (A 7 )a 7 ] . (66) 

In numerous experiments with quantum dots the temperature is usually lower 

than the level spacing. At this condition, if the transition energy A 7 is above both, left 
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and right electrochemical potentials, s 7 is exponentially small in Eq.([66p. This makes iV 7 
exponentially small as well. If, vice versa, A 7 is below both /^ >7 ., s 7 ~ 1, and we obtain 
JVo < 1. Thus, the only case, when N and N p are finite, is when A 7 is located within 
the "conducting window" (CW), i.e., in the interval of energies between \x T and /xj. This 
conclusion is valid, of course, for a finite width T 7 of transitions A 7 , when s 7 is defined by 
Eq.fU). 

Let us introduce the function 

d> 7 = (s- 1 - l)" 1 (67) 
that couples the population numbers of empty and single-electron states: 

iV 7 = iV $ 7 . (68) 

The normalization condition yields 

N ° = ^ = rrfei;- (69) 

In order to display highly non-linear dependence of the population numbers on the applied 
voltage, let us consider the limit s 7 — > s 7 . As a result, we obtain 

$ 7 _ $ 7 = e-P^-ti x (70) 

1 + e ^(A 7 -/,) [ CQsh (/3eV/2) - 5a sinh (/3eV/2)] 
x e /3(A 7 -„) + cosh (/3eV/2) + 5asinh(/?e^/2) ' 

where we have introduced the following notations: /i = (/^+ A t r)/2 , = /x r and the 
degree of asymmetry of the contacts 5a = a 1 — aL. 

At eV = we have $ 7 = exp{— (3 (A 7 — //)} and, therefore, the Gibbs-ensemble limit 
holds: 

1 e -/3(A 7 - M ) 

= i + J] e -P{^-^) ' ^ = 1 + ^ e-^-^' ^ 
In the limit of a small bias voltage, (3eV/2 <C 1, we obtain 



$ 7 « e" 



(72) 



2 

The most prominent phenomenon occurs for a symmetric coupling to the contacts, 5a = 0. 
Namely, in this case we have 

* _ c -ff(A-v-M) 1 + e^-d cosh (PeV/2) 
7 e /?(A-y-„) + cosh (peV/2) ' 1 J 
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When a "level" A 7 gets into the CW, 5eV = {fi r < A 7 — // < //;}, 
exp {f3 (A 7 — + el^/2)} ^> 1. At large bias voltage $ 7 — > 1 and all population 
numbers become equal to each other 

No = N, = — =J (74) 

Population numbers of the states with the energies outside of the interval 5eV are equal to 
zero. Thus, we find a remarkable feature of the correlated transport for a symmetric-coupling 
design of the device: large enough bias voltage equalizes the population of conducting states. 
We will return to this feature later. 



B. The real-time equation for energy shifts. 

We have to define the renormalized transition energy A 7 which is still unknown. The 
transiton-energy shift ()38|) contains the integral 

l 7 = i rftiy 7 (t,ti)D 7 (t 1? t + ), (75) 

J t 

The continued to real time expression, corresponding to Eq. ([75]). will be denoted as Z 7 . We 
have 

K tft = T, L n-^ (76) 



where (see Ref 



7i 



l,= f dr{V>(t,r)D<(T,t)-V<(t,r)D>(r,t)} 

J — oo 

dE ILy(E) 



2mE + i5'' 



(77) 



and 



/OO J 
^ {V>( U )D<( U + E) — V<(u)D>(u + E)} (78) 

Since we assume the adiabatic switching of the interaction and, therefore, ignore the con- 
tribution from the part of the contour (to, to — i0) at to — > — oo, the equivalent forms of 
this expression can be obtained by either consideration of the integral (J75|) on the Keldysh 
contour, or by application of the Langreth' rules (in this case lim#_^+o W-y^i ti)D^(ti, t')] < = 
lim t ^ + o [V<(t,t 1 )I^(t 1 ,t^+V^(t,t 1 )D<(t 1 ,t') 
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Thus, substituting Eqs.(J3SJ) into Eq.(|7Hj). we obtain 



dE ILy(E) 
2iri E + i5 

dE D<(E) 



2jci E — uj + i8 



The GFs D< have a physical meaning of distribution functions. Since we consider wide 
conduction bands in contacts, the width functions V are replaced by a constant and extracted 
from the integral in Eq. (j79j) . In fact, the limits of the integration are determined by the 
bottom and the top of conduction bands ("left" and "right"). It follows from the definition 
of the width, Eq.(|5J): T^(uj) ^ only within the interval of u, where the density of states 
of the corresponding conduction band is nonzero. The remaining integral describes the 
distribution of the density of the mixed charge ((c^Xj) over the energy, which is responsible 
for the renormalization of the transition energy A 7 . 

Inserting into Eq. (|79jl the expressions for "lesser" and "greater" locator, which within 
our approximation are D < = 2iriL^f and D > = — 27uL 7 (l — /), we obtain a simple formula 
for the partial shift of the transition energy 

1 r^.r^r ,„r 7 M[/(i?)-/»] 



L = - duo dELJE) 7V ' ,; v /J (80) 

7 vr J_ w J_ w 7V E-u + %5 

Here, W is a half of bandwidth, which is the same in the both contacts. In fact, the integral 
limits can be replaced by infinities, since the bandwidth is much bigger than any other 
parameter in the theory. 

The right-hand side of Eq. (|79j) contains dressed locators . Therefore, the system of 
equations 

A 7 = A 7 °) + Sf /4 , 

e** = x;^-^ ( 81 ) 

71 

and Eq.(jZHJ) have to be solved self-consistently. 

There are a few remarks in order. First, the Lorentzian width, Eq. (|58jl . is determined by 
the product P 7 T 7 = (N + N^) T 7 . Therefore, if the transition energy A 7 is outside of the 
resonant CW, the population numbers are almost equal to zero. In this case, the coupling 
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between the QD and contacts does not affect the discrete transitions (their widths are zeros). 
Second, the imaginary part of Eq. (JHUj) is exactly equal to zero. Thus, this expression provides 
a simple shift of the transition energy (as it should be). Third, using the definition of the 
width T p , Eq.(jHJ), one can write Eq.fjHUJ) in the form 



1 



oo 



f x (E) - f(e Xkcr ) 



x l--ImD«(E)j , (82) 

which reveals the other aspect of the shift. Namely, estimating the integral at low temper- 
atures, we obtain 

^-^(l</ + K/) (83) 

R = 



[(A 7 - ftif + (r 7 p 7 ) 2 ] [(A 7 - fi r f + (r 7 p 7 ) 2 ] 

Here, we introduced the mixing |vo i7 | 2 = ^2kae\ I v ka,-y 1 2 anc ^ a weighted density of states 
g°(c<j) = J2kae\\ ^ | ^fccr.-y | 2 / l^o^r) fi( u ~~ £ \ka)- Evidently, that the smallness of the coupling 
constant \vq /2W may be compensated by a large logarithm, when the transition energy 
is in the proximity of one of the electrochemical potentials. On the one hand, it may be 
strong and, moreover, it is sensitive to the bias voltage, since /i// r = fi±eV/2. On the other 
hand, the effect is logarithmically weak and, besides, the infra-red cut-off in the integral 
is max{TP, T}, i.e., with an increase of the temperature the renormalization becomes less 
effective. However, the effect is quite appreciable numerically (see below). Thus, we conclude 
that the attachment of contacts to QD shifts the transition energies logarithmically via 
formation of mixed charge (c^Xj in each available channel 7. 



C. Current in the diagonal approximation 

Now we have all ingredients for calculation of the current. The final expression for the 
current is obtained by the substitution Eqs. (j57|) . (f59|) . (}60|) into Eq.((ZJ): 

Ji = ^Ys P -i j dur 7 a l ja ;[fi(uj) - f r (oj)]L 7 (u). (84) 
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This expression is almost identical to the one of the non-interacting electron problem. The 
only difference is that the expression for the SCI regime contains the effective width TP, 
where = N + iV 7 is a combination of population numbers. All P 7 < 1 due to the 
normalization condition, Eq. (J63|) . At Tf3 <C 1 the integral can be estimated by replacing 
the Lorentzian with the delta-function. As a result, we obtain an analogue of the Buttiker's 
formula for the strong Coulomb interaction regime 

T-—X^PT 1 - sinh [eV/2T] 

h ^ 7 cosh [(A 7 -^)/T] + cosh [ e y/2T]' 1 j 

Here, the coefficients a 7 r determine a relative transparency of the "left" /"right" junction. 
The degree of the channel opening is defined by the product r 7 a; 7 a; 7 P 7 . The population of 
each "level" A 7 depends on the populations of other levels due to the sum rule, Eq. fjUHj) . 
Below we apply our theory to describe the magneto-transport through a small vertical 
quantum dot under a perpendicular magnetic field. 

IV. PARABOLIC QUANTUM DOT IN MAGNETIC FIELD 

Shell effects are among the most remarkable phenomena observed in vertical quantum 
dots [12]. In virtue of the potential symmetries, the orbital motion of electrons could lead to 
the degeneracies even at strong intra-dot Coulomb interaction (see, for example, a discussion 



about hidden symmetries for a parabolic potential in Ref 29). Numerous papers devoted to 
the analysis of nonlinear transport through QDs in magnetic field are focused, however, 
on the effect of spin splitting of levels in the dot. Simple estimations evidently indicate 
that the orbital effects are much stronger than spin effects and as yet are not discussed 
in literature related to quantum transport. For example, in a two-dimensional harmonic 



oscillator model for QDs upon the perpendicular magnetic field (cf Ref|30j) the effective 
spin magnetic moment is fi* = Ql^b with fis = \e\h/2m e c and the effective Lande factor 
gi = 0.44. For GaAs the effective mass m* = 0.067m e determines the orbital magnetic 
moment for electrons and gives fig ~ 15/zb, which is 30 times stronger than the spin one. 
In this section we discuss the effects of the orbital motion on the nonlinear (with respect to 
the bias voltage) transport in both the WCI and SCI regimes. 

We consider a dot with a circular shape, u x = u> y = u in a perpendicular magnetic 
field B (Ref. ,30). The dot eigenmodes are Q ± = {Q±uj c /2) with Q = + ^2/4. Here 
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uj c = ^!-B. We choose the position of the QD potential well so that the first level in the 
QD is above /x. In this case at zero bias voltage the QD is empty. 

As was shown above, the coupling to contacts leads to the shift of the transition energies 
and non-Gibbs behavior of the population numbers of QD states. In addition, the spectrum 
of QD in a magnetic field displays degeneracies. Let us illustrate these features in a simple 
example of two degenerate transitions |0) — > \p 2 ). The solution of self-consistent equa- 
tions for this case is displayed on FigO The QD with bare energies is empty (iVo = l)until 
eV riches the bare level e p /uo = 1.5 that is above the Fermi energy without the shift. How- 
ever, the renormalized single-electron states are filled at much smaller values of the bias 
voltage (see solid lines), since the corresponding transition energies (A lQ /u = 0.5) are be- 
low the chemical potential. One observes also that the population numbers are equalized 
(=1/3) after the voltage riches the corresponding transition energy, in the both cases. Below 
we will demonstrate how this mechanism shows up in the transport properties. 
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FIG. 1: (Color online) Population numbers Nq,N p for three states: empty state |0) and two 
degenerate ones \pi), \p2), - as a function of the bias voltage. The energies are given in units of 
confinement energy ojq. The electrochemical potential is chosen at n/uio = 1.0 and T/loq = 0.033. 
The solid (dashed) line displays the behavior of the population number with the renormalized 
(bare) energy. 

A. Shell effects 

The response of the QD to applied bias voltage is determined, evidently, by a electronic 
levels within the conducting window (/i;, /i r ). The number of levels in the window depends : i) 
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on the relative position of the bottom of the confinement potential in QD and the equilibrium 
position of the electrochemical potential (eV = 0); ii) on the strength of magnetic field. 
Although the transport properties in the WCI regime are known, we shortly summarize the 
corresponding effects in order to compare it with the ones in the SCI regime. 

The formula for the current for the WCI regime has the same form as Eq.([84j). but 
contains normal width T instead of effective one TP. Thus, putting formally P = 1 in 
Eq.flBU), we obtain: 



J 



WCI 



e 7 p + [r»r 

In the wide-band mentioned above, the width does not depend on energy, T 1 ^ "(a?) 

T l J, r . Then at small dot-lead coupling we find: 



(86) 



Ji 



4ire sr~^ T^r^ 

i\ 



[fifa) - /r(e 7 )]- 



(87) 
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FIG. 2: (Color online) Magnetic field dependence (t = lo c /ujq) of: a) a Fock-Darwin spectrum; 
the position of the chemical potential fx is displayed by a dashed line; b)the tunnelling current J 
through the quantum dot; c)the conductance G. See text for the definition of J and G. All energies 
are given in units of the confinement energy ojq 

The shell structure in the spectrum of a QD is evidently observed in the WCI regime (see 
FigJ2t)- We use a fixed value for the confinement energy uu . The validity of the assumption 
about the independence of this energy on the number of electrons in the dot, actually, is 
not obvious at all. We believe, however, that strong enough external confining potential and 
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a small number of electrons (i.e., < 20) all results would persist, since the magnetic field 
contributes additionally to the external confining potential. 

At given /i = (//£, + //r)/2 (that simulates in our case the choice of the gate voltage V gate ) 
the number of conducting channels is determined by the number of the Fock-Darwin levels in 
the conducting window < AE < /i£ (CW). This number is changing when the magnetic 
field grows: the levels with higher values of orbital (and spin) momenta m move down faster 
than those with lower momenta and, therefore, the high-lying levels with large values of 
L z = m should unavoidably show up in the CW at large enough magnetic field. The levels 
with the negative orbital momenta L z = —mf move up and leave the CW, decreasing, thus, 
the number of conducting channels. These two processes result in oscillations in the current 
(see FigJ2b)- Thus, each new level entering the CW determined by fi — eV/2 < e 7 < fi+eV/2, 
produces step in the current with the height Jo,7 — ^ireT^T^/ (ftT 7 ). In a narrow energy 
interval the widths T 7 differ from each other only slightly, Jo i7 ~ Jo = weT/h, then the 
"reduced" current J/ Jo should display integer steps. If some level e 7 is n 7 times degenerate, 
the step increases n 7 times. The effect of degeneracy of the spectrum becomes transparent 
at the special values of the magnetic field 

t = u c /u = {k- 1)/Vk, (88) 

which are determined by the ratio = k = 1,2,3,... of the eigen modes of the 

two-dimensional quantum dot (see Fig|3j). 

At these values the current drastically increases, whereas for the values which are slightly 
larger than t , the negative differential conductance, G = dJ/dV, arises. This is illustrated 
in figEfc. 

Fig0] displays the renormalized spectrum, obtained from the self-consistent solution of 
the system of Eqs.(|81 )l .(|82 )l . As expected from the estimation (|8H|) . the renormalization is 
quite large. Due to a large number of the transitions, contributing into the energy shift 
(see Eq. fl82|0 . all transition energies are shifted more or less homogeneously. Thus, in spirit, 
this is, indeed, the MFA: the more transitions contribute into renormalization the closer is 
the shift of each level to an average value (mean field) (see Fig0J). The most prominent 
feature caused by these shifts of the transition energies is a decrease of the threshold of the 
bias voltage to observe the non-zero current. At small eV and chosen fixed \i (or the gate 
voltage) the first transition A 10 < /i r , and the current, Eq.(|51j). is zero, in spite of P 0j ^ 
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FIG. 3: (Color online) Magnetic field (u c /u)q) dependence of: the Fock-Darwin spectrum (in 
units ujq) (left); the tunnelling current Jwci (hi units Jo = T/4/i) through the quantum dot 
(right). Arrows indicate the cuts in the planes that correspond to the fields uj c /^o = 0.0, 0.34, 0.7, 
respectively. The degeneracy at uj c /ujq = 0.0, 0.7 is clearly seen in jumps of the current at T = 0. 




Magnetic field, to c /i>3 Q 

FIG. 4: (Color online) The bare (dashed line) and renormalized (solid line) energies as a function 
of the magnetic field. The parameters of calculations are the same for the WCI and SCI regimes. 
Five lowest level of the Fock-Darwin spectrum were taken as the bare energies. 

(No = 1) (see Fig|SJ). At higher voltages the CW contains n w electron states and, according 
to Eq. (p3§j) . P 7 = 2/(nvi/ + 1). As a result, the SCI current is Jsci = 2Jo n w/( n w + 1), 
where J = ire Y 1 j [hnw] = neT/ (h). The WCI current, however, is Jwci — Jo n w- 
Thus, even for a large bias voltage eV, the SCI current is weaker than the WCI one, by 
factor of r] = Jsci/Jwci — 2/( n w + 1) (until the other charge sector is not switched on at 
eV ~ U). 

Let us compare the manifestation of shell degeneracies denned by the condition, Eq.(|8H|). 
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FIG. 5: (Color online) Population numbers N p (left panel), the current Js/Jo and the conductance 
Gs = dJs /d(eV) as a function of the bias voltage (right panel) in the SCI regime at zero magnetic 
field. In each shell k we have 2(k + l) degenerate orbitals characterized by the same N p (p = k + 1). 
In particular, iVo =< X 00 >, N\ = iVfc = o )7o =i = iVfc = o j7o= 2 etc, where is the orbital index in the 
shell k. Arrows at the bias voltage axis indicate the position of the bare e\ and the renormalized 
Aio energies from the shell k = 0. Their exact values and the parameters of the calculations are 
displayed on the left panel. On the right panel, the rational numbers characterize the height of the 
n-th step in the SCI current, J sa/ Jo = 2(n + l)(n + 2)/(n 2 + 3n + 3), for the last filled shell n. 

in nonlinear transport at the WCI and SCI regimes. We consider first r = 1, i.e., a zero 
magnetic field (see FigEI) - In this case, each shell k has the degeneracy g%. — k + 1. If in 
the transport window the last shell n is filled, the total number of states involved into the 
transport is riw + 1 = 2 X]fc=o(^ + 1) + 1 = ^ 2 + 3n + 3 (2 is due to the spin degeneracy). 
Consequently, the height of the n-th step in the SCI current is Jsa/Jo = 2(n + l)(n + 
2)/(n 2 + 3n + 3), which is smaller than the WCI current by factor rj = 2/(n 2 + 3n + 3). 
Since $ 7 ~ e -£(A 70 -M) [l + a ( e \/) 2 ] at small eV (see Eqs.dZOI)), these effects can not be 
seen in the linear conductance. Another effect (which is not seen in the master-equation 
approach) follows from Eq. (J%Tj) (see also FigJTJ): the coupling pushes the transition energies 
A 7 o down compared to the bare energies e 7 , which decreases the bias voltage threshold, 
as was mentioned above, for the current to be observed. At r = 2 {uj c /uoq ~ 0.7), a new 
shell structure (see FigOJ) arises as if the confining potential would be a deformed harmonic 
oscillator without magnetic field. The number of levels are just the number of levels obtained 
from the 2D oscillator with uj > = 2uj < {uj > and u < denote the larger and smaller value of 
the two frequencies). In this case nw — (n + 2) 2 /2 if the last filled shell is even, and 
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n w = (n + l)(n + 3)/2 if it is odd, and these numbers define the heights of steps in both 
the WCI and SCI regimes. 
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FIG. 6: (Color online) The magnetoconductance G = dJ/(eV) at the WCI (left) and the SCI 
regimes (right), a) The contour plot for the WCI conductance displays the Fock-Darwin spectrum 
of the parabolic QD in the magnetic field. b)The conductance in the WCI regime, at fixed values 
of the bias voltage. The results for eV/wo = 1, 5, 9 are presented by dot-dashed, dashed and solid 
lines, correspondingly. c)The contour plot for the self-consistent conductance at the SCI regime. 
d)The self-consistent conductance at the SCI regime, at fixed values of the bias voltage. The dot- 
dashed, dashed and solid lines present the results for eV/ujQ = 1,1.25,1.45, correspondingly. The 
parameters of the calculations are displayed in the panel " c" . 

These features of the spectrum can be traced by means of the conductance measurements 
(see FigEJ). In particular, the results for the differential conductance dJ/d(eV) at the WC '. 
regime (FigJBk) resembles very much the experimental conductance discussed in Refslll. 
The increase of the bias voltage allows to detect the degeneracy of quantum transitions 
involved in the transport by means of the oscillations in the conductance (see FigEb)- The 
large conductance magnitude at zero magnetic field corresponds to the large number of the 
levels involved in the transport. With the increase of the magnetic field at fixed bias voltage 
(solid line), the conductance oscillates due to appearance and disappearance of the shells. 
As was mentioned above, the increase of the magnetic field brings into the CW states with 
a large magnetic quantum number m = n_ — n + and pushes up the ones with — m. This 
gives oscillations in the current (see FigOD), which are suppressed, however, together with 
the current itself in the SCI regime (see FigEt). 
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Let us discuss once more the physical origin of the current and the conductance suppres- 
sion by the correlations. An increase of the bias- voltage window involves new and new levels 
(or, more accurately, transition energies) into transportation of electrons through QD. In 
SCI regime (contrary to the WCI one) all processes are developing under some constraints. 
First one is that QD cannot contain more than one electron due to a large energy gap 
between single- and two particle states (our model consideration is performed in the limit 
U — ► oo). Second constraint comes from the normalization condition: each channel partici- 
pates in the transport with the spectral weight P 7 = N + iV 7 < 1, whereas the population 
numbers fulfil the sum rule = N + J^ 7 A" 7 = 1. This automatically leads to the conclusion 
that the more levels are involved, the less contribution each channel makes into the current 
and, correspondingly, to the conductance (see FigEJ). Furthermore, since the parabolic QD 
has a spectrum with the shell structure, the number of involved levels grows very fast. This 
explains why a strong suppression of the transport by the correlations occurs in very narrow 
interval of the bias voltage (see FigEJi). Note that in spite of the strong suppression of 
the SCI conductance by the intra-dot correlations, it is possible to observe a fine structure 
produced by the shell structure as well at a certain values of the bias voltage (see FigEJl). 

V. DISCUSSION AND SUMMARY 

In this work we perform an analysis of the non-linear transport through QD, attached 
to two metallic contacts via insulating/vacuum barriers, in the regime of strong intra-dot 
Coulomb interaction. The calculation of the current and conductance is based on the for- 
mula, derived by Meir and Wingreen j^J. It expresses the current in terms of Green func- 
tions of fermions in a QD. In the regime of strong Coulomb interaction the fermion variables 
are not good starting point for a constructing the perturbation theory with respect to the 
small transparency of the barriers. For this reason we use the approach where an active 
element (quantum dot) of the device is supposed to be treated as good as possible (ide- 
ally, exactly), whereas the coupling to the metallic contacts is treated perturbatively. This 
leads to description in terms of many-electron operators; particulaly, here we find conve- 
nient to use the Hubbard operators. We have re-written the Meir-Wingreen's formula in 
terms of Green functions (GFs), defined on Hubbard operators. The Wick's theorem does 
not work for such GFs. In our approach this problem is solved by means of the extension 
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of the diagram technique for Hubbard-operator Green functions (developed earlier for a 
equilibrium in Reflisl) to the non-equilibrium states. In the first step, the exact functional 
equations for the many-electron GFs are derived and used for the iterations in the spirit of 

n 

the Schwinger-Kadanoff-Baym approach (see also Refll8l). This gives rise to the diagram- 
matic expansion. In the second step the analytical continuation of these equations to the 
real-time axis transforms the system to the non-equilibrium form. 

At the limit of infinitely strong intra-dot Coulomb interaction, the lowest approximation 
with respect to a small transparency of barriers happen to be the mean field theory. In 
this theory the population numbers of the intra-dot states and the energies of the transition 
between these states have to be found self-consistently in non-equilibrium (but stationary) 
states. Thus, the non-equilibrium mean field theory for multiple-level QD has been con- 
structed for non-zero temperature (to the best of our knowledge this is the first attempt in 
literature). It worth to note that a common practice is to consider the population numbers 
of many-electron states in the Gibbs' form (following, e.g., RefJ^). In such an approach the 
population numbers do de pen d on the magnetic field, however the dependence on the bias 
voltage is missing 

(cf RefsflHS- 

One of the questions we have addressed in the paper is to what extent the levels (energies 
of the transitions) of the QD are sensitive to the applied bias voltage eV. We have found that 
this dependence is logarithmically weak (see Eq l83|) . In contrast, the population numbers 
manifest much stronger and highly non-trivial dependence on the bias voltage. At zero 
bias voltage they display the standard Gibbs' statistics. However, at the voltage higher 
than a certain critical value (determined by the position of renormalized transition energy) 
the population numbers start equalizing! This occurs with only those of transitions which 
have energies within the conducting window (fi r ,fii). This result is obtained analytically 
in simplified equations and confirmed by the numerical solution of a complete system of 
non-linear integral equations. 

The expression for the current for SCI regime that we have derived (Eq. (j84jl ). demon- 
strates that the degree of the nonlinearity of the transport is essentially determined by the 
spectral weights P p of the dot states. Since they are just the sum of the population numbers, 
Pp — N + N p , the non-trivial behavior of the latters (see Eqs.(|M jl .(|H7 jl .(|Hn )l ) is reflected in 
transport properties via a spectral weight of each conducting channel. 

Applying further these equations to the problem of the transport through a two- 
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dimensional parabolic quantum dot, we obtain a simple expression which relates the height 
of the n-th step in the current to the number of states participating in the transport. In the 
WCI regime, at specific values of the magnetic filed, Eq.()88|). we predict a drastic increase 
of the current through the dot due to the shell effects. Possibly, this effect can be used 
in applications (i.g., for magnetic-field sensitive switch on-off devices). In the SCI regime 
the transport is strongly suppressed. The origin of it again is connected with the spectral 
weights, regulating the contribution of each conducting channel. First, the sum of popula- 
tion numbers is equal to one. Second, the number of levels within the conducting window 
grows very fast due to existence of the shell structure in the spectrum of parabolic QD. 
Third, in the limit U — > oo the QD cannot contain more than one electron. Therefore, the 
contribution of each channel unavoidably drops down when expanding by increasing voltage 
conducting window embraces more and more transition energies. As a result, the current 
tends to a certain limit, whereas the heights of the peaks in the conductance tend to zero. 
Thus, we conclude: while in the WCI regime the fine structures in the conductance within 
one "diamond" can arise due to re-scattering/co-tunneling processes, in the SCI regime it 
has the other origin. Namely, the role of single-electron levels is played by the energies of 
the transition between many-electron states, and the contribution of each transition into 
transport is regulated by the strength of its coupling to the conduction-electron subsystem 
and the spectral weights. 
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APPENDIX A: RELATIONS BETWEEN FERMI AND HUBBARD OPERA- 
TORS 

It is convenient to use different notations for Hubbard operators that describe Bose-like 
transitions (without change of particle numbers) and Fermi-like ones (between the states, 
differing by one electron). The former and the later ones are denoted as (Z) and (X), 
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respectively: 

= H (1 - = |0> <0| , (A.l) 

71 

ZTi= n (l-< 77l )4d 7 = | 7 )( 7 |, (A.2) 
X° 7 = n (1-^)^=10X71, (A.3) 

71^7 

Z 00 |0) = |0), \j) = \j) , (A.4) 

where Eq. (|A.l|) defines the vacuum state and the product is taken over all single-particle 
states 71 of the QD confining potential. In this case the expansion of unity acquires the 
form: 

i = z m + J2 zl1 - ( A - 5 ) 



7 

The expectation values N = (Z 00 ) and iV 7 = (Z 77 ) with respect to state/ensemble of 
interest are population numbers of corresponding many-electron states. The Hamiltonians 
of closed QD in these terms become: 

Hqd 1 = E £ A d ^ H QD = + £ e y Z^. (A.6) 

7 7 

The WCI regime corresponds to U — > 0, while the SCI regime takes place at U — > oo. The 
energy of empty QD £ is determined by the potential depth and can be chosen arbitrarily. 
While the eigenvalues e 1 are defined by shape of the potential, the transition energies are 
= e 7 — eq. We recall that the position of the single-electron transitions with respect to 
the electrochemical potentials is essential for the transport properties of the system. 

To write H t in these variables, we have to express the annihilation (creation) operator 
d 7 (dl) in terms of Hubbard operators: 

d 7 = i • d 7 ■ 1 = (0| d 7 | 7 ) X° 7 . (A.7) 

Of course, the expansion (jA.7|) should contain the transitions from one- to two-particle, 
from two- to three-particle and so on configuration A n 

d 1= (A (n) ,n|rf 7 |A(" +1 ),n + l)X A( " )A(n+1) , (A.8) 

{A,n} 

that defines the anticommutator relation 

i = K.<}= }2 zAMA(n) - (a. 9) 

{A,n} 
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However, in the limit U — > oo the above equation is reduced to Eq. (|A.5|) . In other words, 
we consider the transport within the first conducting "diamond" only, in coordinates of 
gate and bias voltages. The H t remains unchanged in the WCI regime, whereas in the SCI 
regime it can be written in the form 

H t SCI = E (<A^° 7 + H {A.10) 

fc,frgA,A,7 

where 

<r; 7 = £L; 7 (0Kl7>, 

v^ M = v^\dS r \0) (A.ll) 

The matrix element v ka . Q ^ contains information on specific features of attachment between 
a quantum dot and contacts. Here, X = l,r denotes the left and the right lead, respectively. 
One may notice that Hubbard operators take into account corresponding kinematic restric- 
tions placed by the strong Coulomb interaction. The prize for this is, however, non-trivial 
commutation relations: 

{X°\X-<'°} = 5 y yZ 00 + Z"i'\ [Z 00 ,X^] = X°\ 

y z ii> ; x o 7 ] = x o 7 " _ (A. 12) 

Thus, our model Hamiltonian takes the form 

H = EkaeX,X £ ^ac[ ka C Xk a + B Z 00 + J^j^ 2 ^ 

+ E fcCTe A,A, 7 (< 7 4^° 7 + h.c) (A. 13) 

We assume also that {c X ka, c[, k , a ,} = 6 k> k>8a,*'$\,\', {c\ka, d\} = 0. 

As it was discussed in the Introduction, the particular sector of the transitions (0 <^=^> 1) 
is remarkable due to the fact that the bare energies e 7 and the tunnelling matrix elements 
coincide in both the WCI and the SCI regimes. This opens a nice opportunity to study the 
role played by the strong Coulomb interaction in formation of transport properties in a most 
refined form. Indeed, in the SCI case the electrons in QD can "see" the conduction electrons 
in the attached contacts only through the prism of the in-dot kinematic constraints. This 
results in the renormalization of the transition energies. 
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APPENDIX B: CURRENT 



Here we re-derive the expression for the current obtained by Meir and Wingreen 23j in 
terms of Green functions (GFs) for the Hubbard operators. Thus, 



E {<,(4 k ^)-v\ M (X^c lk ,)} 

E / ^{<,7 G 7V^)-<^,7M} 



ie 



e 
h 



(B.l) 



where the matrix elements v l ka (v, 1 ka ) are defined by Eq. ()A.ll|) . Equations for the Green 
functions G^^a and Gk a ,j on the Keldysh contour are 



G iM (t, t') = j c dnG^y (r, n V, , (^(ri, t') 



G fcCT , 7 (r, r') = / dnCf} (r, n)^ CT ,Gy i7 (r, n) 



(B.2) 
(B.3) 



Here C^ a (ti, r') is a bare GF of conduction electrons in the left lead. Expressing GFs G via 
integrals on real-time axes and making a Fourier transformation with respect to a difference 
of times (the steady state is considered only!), we obtain: 

G< k M = G< 1 ,{u J )v\ IM C^ A {u) + 
+ cT{u)v>G<(u) 



Here, the Fourier transformation is defined as: 

duo 



G(t)= I ^e- i " t G(u),G(uj)= I dte luJt G(t). 



(B.4) 



(B.5) 



Inserting Eqs. ()B.4jl into Eq.(jBJJ), we have: 



rl,</ 



J du v£{u)G*Au) + V^ R {u)G<{u) 



(B.6) 
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Here, we have introduced the effective interaction with conduction electrons in the lead 
A = /: 

kaeX 

where r\ = {R, A, >,<}. Since in the steady state the current is homogeneous, i.e., J; = 
(J/ + Ji)/2 = (J; — J r )/2, we can symmetrize Eq. (jB.6|) : 

Ji= ^ / <M(^» - v#») - 

-W?>) - ^>))]G< y (u;) + (B.8) 
+^ / dw(^<( w ) - ^/») (G£,(u;) - G^( W )) 

Let us calculate the effective interaction VlfJu). The conduction-electron "lesser" and 
"greater" GFs for the left lead (A = /) or for the right one (A = r) are: 

CS K (M') = i(c3^(Oc>te(t)) (0) = <e^^)/ A ( eAfc ) 

cJfftO = -<(ca*.(0cUo) <0) = 

Cf K M = • | d(t-f)e i(w - eA *- )( *^ ) /A(eAto) = 

= 2ttz5 (a; - e\ ka ) f x {uj) 
C£>(a;) = -i y d(t - Oc'^^^AC-eAto) = 

= -2tu5 (a; - £ A fe CT ) f\(-w) 

For the retarded and advanced GFs we have 

= ~iO(t-lf)({c^(t),c]^(l/)})^ = 
= -i6(t - £') e -fe*U*-*') 

C£ )fl (w) = -* / d (* - 0^(* - tV^^^' = 
1 

c^ A (t,t>) = tf(-t+o({cA*(*),4,(0})f = 

C<° )j4 (u;) = ^ y d(t-t')0(-t + t') e ^- e ^)(*-*') = 
1 

u - £\ka - iS' 
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As a result, we have 



and 



X,R/ 

W)J — 



kcreX 



fc(T,7 



kaeX 



With the aid of above equations we can rewrite Eq.(B7) in terms of the width functions: 
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